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1. Introduction 

The field of Numerical Relativity (NR) has progressed at a remarkable pace since the 
breakthroughs of 2005 P, El E] with the first successful fully non-linear dynamical 
numerical simulation of the inspiral, merger, and ringdown of an orbiting black-hole 
binary (BHB) system. In particular, the 'moving-punctures' approach, developed 
independently by the NR groups at NASA/GSFC and at RIT, has now become the 
most widely used method in the field and was successfully applied to evolve generic 
BHBs. This approach regularizes a singular term in space-time metric and allows the 
black holes (BHs) to move across the computational domain. Previous methods used 
special coordinate conditions that kept the BHs fixed in space, which introduced severe 
coordinate distortions that caused orbiting-BHB simulations to crash. Recently, the 
generalized harmonic approach method, first developed by Pretorius pQ, has also been 
successfully applied to accurately evolve generic BHBs for tens of orbits with the use of 
pseudospectral codes [U [3] . 

Since then, BHB physics has rapidly matured into a critical tool for gravitational 
wave (GW) data analysis and astrophysics. Recent developments include: studies of the 
orbital dynamics of spinning BHBs [HI EJ EJ [101 EEU 12] , calculations of recoil velocities 
from the merger of unequal mass BHBs [T3l [HJ US] , and the surprising discovery that 
very large recoils can be acquired by the remnant of the merger of two spinning BHs 
[161 H3 UHl H3 [201 [2U [221 [231 I2H [2S1 [261 [2ZI 1211 El [291 EQl [31], empirical models 
relating the final mass and spin of the remnant with the spins of the individual BHs 
[32] 1331 IM1 1351 136] 1371 EH E2] , comparisons of waveforms and orbital dynamics of BHB 
inspirals with post-Newtonian (PN) predictions 0QHH1H21S31IS1H51HS1S7J, and 
simulations reaching mass ratios [IE] q = 1/100. 

1.1. Kicks and Super-kicks: A brief history 

The first in-depth modeling of the recoil from the merger of non-spinning asymmetric 
BHBs was done in Ref. [15], where it was shown that the maximum recoil is limited to 
pd 175 kms -1 . Soon after, other groups showed that the maximum recoil for spinning 
binaries is much larger. In Ref. [16] and [30] (which were released within a few days of 
each other), it was shown that the maximum recoil for a spinning binary with one BH 
spin aligned with the orbital angular momentum and other anti-aligned id 475 kms -1 . 
Within a day of the initial release of the preprint for Ref. [16J, our group released a 
preprint p2] for more generic spin and mass configurations, showing that the recoil is 
considerably larger if the spins are anti-aligned and pointing in the orbital plane. In [17] 
we measured recoil velocities beyond the maximum predicted for the configurations 
in [HUE)]. 

An initial analysis of the results in [T7] indicated that the maximum velocity 
exceeded 1300 kms -1 for spins lying in the orbital plane. This initial estimate was 
based on two 'generic' configurations. A subsequent analysis that incorporated the 
angular dependence of the projection of the spin on the orbital plane, showed that 
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the maximum recoils is ~ 4000 kms -1 jl9]. Based on our conclusions in the preprint 
of [17], the group at Jena performed a set of two simulations in this 'maximum- kick' 
configuration and measured recoils of around 2600 kms -1 . Our modeling showed that 
this recoil velocity actually depends sinusoidally on the angle that the spins make with 
the infall direction. By evolving a set of configurations with spins at different initial 
angles in the orbital plane, we found that the recoil reaches a maximum of ~ 4000 km s -1 
|18| 155] . In Ref. [H5], as part of our analysis, we proposed an empirical formula for the 
recoil based on post-Newtonian expressions for the radiated linear momentum. The 
formula correctly predicts the sinusoidal dependence of the 'maximum-kick' recoil. 

In Ref. [2S] the recoil for unequal-mass, spinning binaries, with total spin equal to 
zero and individual spins lying in the initial orbital plane, was measured. These S = 
configurations are preserved under numerical evolution and lead to minimal precession 
of the orbital pane. Interestingly, they found that the recoil scales with the cube of the 
mass ratio q, rather than the expected q 2 seen in post-Newtonian expressions for the 
recoil. In a subsequent study [19], our group found that the recoil scales as q 2 for the 
more astrophysically important precessing case. 

1.2. Phenomenological modeling of the recoil 

In [T7] we introduced an empirical formula for the recoil, augmented in [39] which has 
the form. 



where rj = q/(l + q) 2 , the index _L and || refer to perpendicular and parallel to the 
orbital angular momentum respectively, ii, e 2 are orthogonal unit vectors in the orbital 
plane, and £ measures the angle between the unequal mass and spin contribution to 
the recoil velocity in the orbital plane. The constants H$ and K$ can be determined 
from new generic BHB simulations as the data become available. The angles, Ga and 
&Si are the angles between the in-plane component of A = M(S2/rri2 — Si/m{) or 
S — Si + S 2 and the infall direction at merger. Phases 6 and G>i depend on the 
initial separation of the holes for quasicircular orbits. A crucial observation is that 
the dominant contribution to the recoil is generated near the time of formation of the 
common horizon of the merging BHs (See, for instance Fig. 6 in ]50j). The formula 
Q above describing the recoil applies at this moment (or averaged coefficients around 
this maximum generation of recoil), and has proven to represent the distribution of 



Vrecoii(<?,a) =v m e 1 + v_l(cos£ ei + sin£e 2 ) + v\\ h\ 




(ct| + q 2 a\) , 



(1) 
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velocities with sufficient accuracy for astrophysical applications. Recently, in Ref. [51] 
a very similar formula was proposed. 

The most recent published estimates for the above parameters can be found in 
[19] and references therein. The current best estimates are: A = 1.2 x 10 4 kms -1 , 
B = -0.93, H = (6.9 ± 0.5) x 10 3 kms" 1 , K = (6.0 ± 0.1) x 10 4 kms" 1 , and £ ~ 145°. 

Note that v» is maximized when q — 1 (i.e. equal masses), the in-plane spins are 
maximal (cti = ot-x — 1), and the two spins are anti-aligned. 

In this paper, we use a new set of simulations of precessing binaries to test the 
recoil formula and obtain a new estimate for K$- 

2. Numerical Techniques 

To compute the numerical initial data, we use the puncture approach |52j along with 
the TwoPunctures [53j thorn. In this approach the 3-metric on the initial slice has 
the form ^ ab = (ipBL + w) 4 <5 a fe, where i/jbl is the Brill-Lindquist conformal factor, 5 ab 
is the Euclidean metric, and u is (at least) C 2 on the punctures. The Brill-Lindquist 
conformal factor is given by ipBL = 1 + Y2i=i "^f /(2|r — fi\), where n is the total number 
of 'punctures', is the mass parameter of puncture i (mf is not the horizon mass 
associated with puncture i), and r*j is the coordinate location of puncture i. We evolve 
these BHB data-sets using the LazEv [54J implementation of the moving puncture 
approach [2j [3] with the conformal factor W = ^/x = exp(— 20) suggested by [H] For 
the runs presented here we use centered, eighth-order finite differencing in space [55] 
and an RK4 time integrator (note that we do not upwind the advection terms). 

We use the Carpet [56] mesh refinement driver to provide a 'moving boxes' style 
mesh refinement. In this approach refined grids of fixed size are arranged about the 
coordinate centers of both holes. The Carpet code then moves these fine grids about 
the computational domain by following the trajectories of the two BHs. 

We use AHFinderDirect |57J to locate apparent horizons. We measure the 
magnitude of the horizon spin using the Isolated Horizon algorithm detailed in |58j. 
This algorithm is based on finding an approximate rotational Killing vector (i.e. an 
approximate rotational symmetry) on the horizon (p a . Given this approximate Killing 
vector (p a , the spin magnitude is 

%] = [ (<p a R b K ab )d 2 V } (2) 

S7r J AH 

where K ab is the extrinsic curvature of the 3D-slice, d 2 V is the natural volume element 
intrinsic to the horizon, and R a is the outward pointing unit vector normal to the horizon 
on the 3D-slice. We measure the direction of the spin by finding the coordinate line 
joining the poles of this Killing vector field using the technique introduced in [8] . Our 
algorithm for finding the poles of the Killing vector field has an accuracy of ~ 2° (see [8] 
for details). Note that once we have the horizon spin, we can calculate the horizon mass 
via the Christodoulou formula m H = ■sj mf rv + S 2 / (Am 2 ^), where m irr = \J A/ (IQti) and 
A is the surface area of the horizon. We measure radiated energy, linear momentum, 
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and angular momentum, in terms of using the formulae provided in Refs. j59j [60] . 
However, rather than using the full ip^, we decompose it into i and m modes and 
solve for the radiated linear momentum, dropping terms with i > 5. The formulae in 
Refs. [59j [60] are valid at r = oo. We obtain accurate values for these quantities by 
solving for them on spheres of finite radius (typically r/M = 50, 60, • • • , 100), fitting the 
results to a polynomial dependence in I — 1/r, and extrapolating to I — [21 EH EH [62]. 
Each quantity Q has the radial dependence Q = Q + IQ\ + Oil 2 ), where Q is the 
asymptotic value (the 0(1) error arises from the 0(1) error in rtp^). We perform both 
linear and quadratic fits of Q versus /, and take Qq from the quadratic fit as the final 
value with the differences between the linear and extrapolated Qq as a measure of the 
error in the extrapolations. We found that extrapolating the waveform itself to r = oo 
introduced phase errors due to uncertainties in the areal radius of the observers, as well 
as numerical noise. Thus when comparing Perturbative to numerical waveforms, we use 
the waveform extracted at r = 100M. 

In order to model the recoil as a function of the orientation and magnitudes of the 
spins, we use the techniques introduced in [19] to locate the approximate orbital plane 
at merger and 3D rotation such that infall directions are the same for each simulation. 
Briefly, this technique uses three points on the trajectories, given by fiducial choices of 
the BH separations, to define the orbital plane and preferred orientation. 

3. Simulations 



Config 


xi/M 


x 2 /M 


P/M 


p 

m\ 




S x /M 2 


Sy/M 2 


Q33TH000 


4.882446 


-1.607923 


0.101163 


0.171173 


0.723529 





0.045983 


RQ33TH000 


4.885558 


-1.609045 


0.101153 


0.171170 


0.723523 





0.045982 


Q50TH000 


4.360493 


-2.163334 


0.119252 


0.230648 


0.618813 





0.082098 



Table 1. Initial data parameters for the quasi-circular configurations. The punctures 
are located at fi = (afi,0, 0) and f*2 = (x2,0, 0), with momenta P = ±(0, P, 0), spins 
Si = (S x , S y , 0), S2 = —qSi, mass parameters m v . The configuration are denoted by 
QXXXTHYYY where XXX gives the mass ratio (0.33, 0.50) and YYY gives the angle 
in degrees between the initial spin direction and the y-axis. In all cases the initial 
orbital period is Muj = 0.05 and the spin of the smaller BH is a = 0.72. Initial data 
parameters for the Q50TH000 and Q33TH000 configurations are given. The remaining 
configurations are obtained by rotating the spin directions, keeping all other parameters 
the same. For the RQ33THxxx configurations, §2 is rotated by 90° with respect to 
the corresponding Q33THxxx configuration. 

We evolve a set of configuration that initially have A = 0, as well as a set of 
configuration with one spin rotated by 90° for mass ratios q — 1/3 and q = 1/2. 
The initial data parameters are summarized in Table [3j Note that the A = 
configurations are unstable in the sense that the system quickly evolved towards a 
nontrivial A. In particular, at merger, where most of the recoil asymmetry is generated, 
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the binary reaches a generic configuration regarding the spin orientations given the 
strong differential precession of each BH spins. 

4. Results 

In Table [2] we summarize the results of the simulations. The table shows the radiated 
energy and recoil (prior to any rotation). Note that these results can be used in 
additional fits of the final remnant BH formulae for the mass, spin, and recoil velocity 
of the remnant, as was done in Ref. [39] using the then currently available results in the 
literature. 



Config 


100(6E/M) 


v x 


Vy 


v z 


Q50TH000 


2.858 ±0.018 


76.89 ± 2.47 


234.40 ±3.12 


167.95 ±3.29 


Q50TH045 


2.817 ±0.017 


-40.47 ± 131.03 


10.46 ± 99.56 


-228.49 ± 5.76 


Q50TH090 


2.778 ±0.016 


189.23 ±0.77 


97.60 ± 0.34 


-543.78 ± 2.29 


Q50TH130 


2.795 ±0.017 


306.54 ± 1.23 


319.10 ±0.02 


-762.32 ± 0.29 


Q50TH180 


2.858 ±0.018 


78.03 ±2.66 


235.61 ±3.24 


-167.47 ±3.37 


Q50TH210 


2.831 ±0.017 


43.42 ± 1.32 


108.60 ±2.94 


131.29 ±3.42 


Q50TH315 


2.831 ±0.017 


295.55 ± 1.26 


341.01 ±0.78 


741.05 ±0.66 


Q33TH000 


1.90243 ± 0.012 


119.18 ±0.36 


212.92 ±2.78 


144.191.08 


Q33TH045 


1.884 ±0.012 


67.29 ±0.59 


132.08 ±2.51 


-134.74 ±0.94 


Q33TH090 


1.878 ±0.012 


02.92 ± 1.87 


132.08 ± 1.51 


-337.04 ± 0.88 


Q33TH130 


1.897 ±0.012 


210.61 ± 1.34 


235.72 ±0.71 


-497.12 ± 1.34 


Q33TH180 


1.902 ±0.012 


119.21 ±0.36 


212.83 ±2.78 


-144.04 ± 1.08 


Q33TH210 


1.889 ±0.012 


73.21 ±0.11 


150.46 ± 2.73 


59.42 ± 1.02 


Q33TH260 


1.878 ±0.012 


73.45 ± 1.79 


127.45 ± 1.58 


272.87 ± 1.14 


Q33TH315 


1.901 ±0.013 


210.14 ± 1.45 


249.05 ±0.87 


489.34 ± 0.76 


RQ33TH000 


1.885 ±0.013 


229.17 ± 1.07 


132.08 ± 1.81 


539.21 ± 1.66 


RQ33TH090 


1.865 ±0.012 


45.62 ±0.11 


162.27 ± 1.00 


-74.423 ± 0.08 


RQ33TH130 


1.860 ±0.012 


106.45 ± 1.16 


82.27 ±0.45 


-427.92 ± 1.43 


RQ33TH210 


1.886 ±0.012 


182.61 ±0.67 


188.75 ±2.00 


-361.19±1.04 


RQ33TH315 


1.861 ±0.012 


121.99 ± 1.42 


79.03 ± 0.45 


462.96 ± 1.68 



Table 2. The radiated energy and recoil velocities for each configuration. Note that 
some of the error estimates, which are based on the differences between a linear and 
quadratic extrapolation in I = 1/r of the observer location, are very small. This 
indicates that the differences between the extrapolation can underestimate the true 
error. All quantities are given in the coordinate system used by the code (i.e. the 
untransformed system). 

Table [3] gives the components of the radiated angular momentum in the original 
x,y,z frame (that of the initial data) using the Cartesian decomposition as in Ref. [60J. 
In Table [4] we give the recoil velocities in a frame rotated such that the orbital 
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Config 

Q50TH000 

Q50TH045 

Q50TH090 

Q50TH130 

Q50TH180 

Q50TH210 

Q50TH315 



SJx 

0.02454 ± 0.00056 
-0.0277 ±0.0041 
-0.05995 ±0.00019 
-0.06155 ±0.00029 
-0.02454 ± 0.00056 
0.00857 ±0.00035 
0.05981 ± 0.00026 



6 Jy 

0.05954 ±0.00012 

0.0572 ± 0.0029 

0.02354 ± 0.00079 

-0.02016 ±0.00047 

-0.05954 ±0.00012 

-0.06326 ± 0.00059 

0.02545 ± 0.00046 



SJ Z 

0.18040 ±0.00200 
0.1808 ±0.0017 
0.1808 ±0.0015 

0.1801 ±0.00164 
0.1804 ±0.0020 
0.1805 ±0.0018 
0.1802 ±0.0017 



Q33TH000 
Q33TH045 
Q33TH090 
Q33TH130 
Q33TH180 
Q33TH210 
Q33TH260 
Q33TH315 



0.0180 ±0.00011 
-0.00930 ± 0.00008 
-0.03151 ±0.00002 
-0.03524 ±0.00047 
-0.01801 ±0.00011 
-0.00013 ± 0.00002 
0.02797 ± 0.00003 
0.03462 ± 0.00041 



0.03134 ±0.00022 
0.03504 ±0.00010 
0.01835 ±0.00013 
-0.00678 ± 0.00002 
-0.03134 ±0.00022 
-0.03608 ± 0.00022 
-0.02366 ±0.00014 
0.00994 ± 0.00002 



0.1309 ±0.0012 
0.1312 ±0.0012 

0.13262 ±0.00081 
0.1325 ±0.0010 
0.1309 ±0.0012 
0.1309 ±0.0012 

0.13254 ±0.00084 
0.1325 ±0.0011 



RQ33TH000 
RQ33TH090 
RQ33TH130 
RQ33TH210 
RQ33TH315 



0.021251 ±0.000001 
0.00462 ± 0.00017 
-0.00981 ±0.00001 
-0.02076 ± 0.00014 
0.01143 ± 0.00001 



-0.00413 ± 0.00003 
0.02225 ± 0.00001 
0.01968 ±0.00017 
-0.00734 ± 0.00008 
-0.01857 ±0.00022 



0.1348 ±0.0012 
0.13479 ±0.00091 
0.1346 ±0.0010 
0.1345 ±0.0012 
0.1346 ±0.0011 



Table 3. The radiated angular momentum. Note that some of the error estimates, 
which are based on the differences between a linear and quadratic extrapolation in 
I = 1/r of the observer location, are very small. This indicates that the differences 
between the extrapolation can underestimate the true error. All quantities are given 
in the coordinate system used by the code (i.e. the untransformcd system). 



plane coincides with the xy plane and the infall direction in this plane is fixed. In order 
to fit the data, we chose K = 5.9 x 10 4 based on previous work [T7J [191 E3] and then 
fit the Q50 simulations for K s , 6 , and 6i. We found K s = -4.2 ± 1.8. While, K 
and Ks are fixed, 0o and 0i depend on the configuration. To obtain these angles we 
fit the the recoil formula for the Q33Txxx and RQ33Txxx configurations separately. 
We then compared the predicted recoil for each configuration with the measured recoil. 
The results are summarized in Table |4j Note that the errors are typically less than 
6 kms" 1 . The largest relative errors are ~ 3.3%, with most errors lying between ~ 1% 
and ~ 2%, rendering the empirical formula accurate for most astrophysical applications. 
We plot the measured out-of-plane recoil and prediction in Fig. [TJ An attempt to fit 
the in-plane recoil produced errors of the order of 50km s -1 . This can be traced to three 
main sources the errors. The uncertainty in how the the recoils produced by unequal 
masses and out-of plane spins contribute to the total in-plane recoil, the error in the 
measurement in the out-of-plane spine due to errors in measuring the orientation of the 
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orbital plane at merger, and, as pointed out in [19] . effects of precession on the in-plane 
recoil. We note that the dominant out-of-plane recoil is well modeled. 



Config 






V 

v z 


V- ( Dredict 1 


error 


O50TH000 


17.47 


-164.42 


-248.44 


-245 851 


2.591 


O50TH045 


-71.85 


-100.95 


196.47 


190.25 


-6.222 


O50TH090 


44.52 


-180.82 


553.49 


551 759 


-1.73 


O50TH130 


31.47 


-242.82 


846.74 


848 283 


1.540 


O50TH180 


18.34 


-165.91 


248.56 


245 987 


-2.575 


O50TH210 


26.93 


-153.45 


-81.52 


-84.3437 


-2.823 


Q50TH315 


28.58 


-242.02 


-832.71 


-834.713 


-1.999 


O33TH000 


117.61 


-157.99 


-203.81 


-208 612 


-3.311 


O33TH045 


110.76 


-132.11 


102.02 


105.811 


1.937 


O33TH090 


122.60 


-120.81 


334.67 


337 689 


-0.893 


O33TH130 


144.38 


-155.39 


549.60 


553 834 


-2.540 


Q33TH180 


117.60 


-157.99 


203.63 


208.286 


3.140 


Q33TH210 


109.72 


-138.43 


-18.17 


-18.7138 


0.501 


Q33TH260 


118.62 


-122.29 


-258.97 


-268.242 


-6.040 


Q33TH315 


144.43 


-160.93 


-546.70 


-551.418 


2.112 


RQ33TH000 


168.76 


-118.44 


-564.09 


-569.559 


-5.469 


RQ33TH090 


116.46 


-133.88 


49.68 


49.8794 


0.203 


RQ33TH130 


121.09 


-92.33 


421.93 


413.897 


-8.033 


RQ33TH210 


155.20 


-147.31 


391.98 


386.733 


-5.246 


RQ33TH315 


125.67 


-89.20 


-460.12 


-464.877 


-4.757 



Table 4. The recoil velocities after rotation and a comparison of the out-of-plane 
component of the recoil (in coordinates aligned with the orbital plane at merger) 
with the predictions of the empirical formula with coefficients K = 5.9 x 10 4 and 
K s = -4.25401. 



5. Conclusion 

We report here on a new set of generic simulations, with no symmetries and different 
choices of the binary's mass ratio and spin direction and magnitudes. We compute 
the radiation emitted by this binaries, and in particular focus on the magnitude and 
direction of the final recoil of the merged BH. While this set of runs can be used for 
fitting to additional subleading terms in the empirical formulae for the remnant BH 
recoil, we choose to use them to test the previously fitted values, extended to include 
effects linear in the spin, as accurate determinations of these parameters will require 
many more simulations. In this paper we showed that the empirical recoil formula 
provides accurate predictions for the recoil velocity from BHB mergers for the new sets of 
configurations. These configurations have fewer symmetries than previous comparisons 
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Figure 1. The predicted and measured out-of-plane recoil for each configuration in 
Table [3] The configurations are labeled 1-20 starting from Q50TH000 as ordered in 
Table H 



and can be considered of generic nature regarding spin orientations and intermediate 
mass ratios and spin magnitudes. This shows that our empirical formula ([!]) can be 
used as a first approximation for astrophysical studies of statistical nature, as we did in 
Ref. [M] and also used for realistic recoil magnitudes and direction when modeling the 
observational effects of recoiling BHs in a gaseous environment such as accretion disks 
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